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ABSTRACT 

We investigate the conditions for the existence of an expanding virial shock in the gas 
faUing within a spherical dark-matter halo. The shock relies on pressure support by 
the shock-heated gas behind it. When the radiative cooling is efficient compared to 
the infall rate the post-shock gas becomes unstable; it collapses inwards and cannot 
support the shock. We find for a monoatomic gas that the shock is stable when the 
post-shock pressure and density obey jcS = {d In P / dt)/{d In p/dt) > 10/7. When 
expressed in terms of the pre-shock gas properties at radius r it reads prA(T)/u^ < 
.0126, where p is the gas density, u is the infall velocity and A(T) is the cooling function, 
with the post-shock temperature T (x v? . This result is confirmed by hydrodynamical 
simulations, using an accurate spheri-symmetric Lagrangian code. When the stability 
analysis is applied in cosmology, we find that a virial shock does not develop in most 
haloes that form before z 2, and it never forms in haloes less massive than a few 
1O^^M0. In such haloes, the infalling gas is not heated to the virial temperature until 
it hits the disc, thus avoiding the cooling-dominated quasi-static contraction phase. 
The direct collapse of the cold gas into the disc should have nontrivial effects on the 
star-formation rate and on outflows. The soft X-ray produced by the shock-heated 
gas in the disc is expected to ionize the dense disc environment, and the subsequent 
recombination would result in a high flux of La emission. This may explain both 
the puzzling low flux of soft X-ray background and the La emitters observed at high 
redshift. 

Key words: cooling flows — dark matter — galaxies: formation — galaxies: ISM — 
hydrodynamics — shock waves 



1 INTRODUCTION 

The standard lore in the idealized picture of galaxy forma- 
tion by spherical infall of gas inside dark-matter haloes is 
that the gas is first heated to the halo virial temperature 
behind an expanding virial shock. It is then supported by 
pressure in a quasi-static equilibrium while it is cooling ra- 
diatively and is slowly contracting to a disc where it can 
eventually form stars. The cooling process thus determines 
important galaxy properties such as the star-formation rate 
and the metal enrichment, so it is necessarily an important 
ingredient in the galaxy formation process. 

However, it is not at all clear that a stable shock can 
persist in the halo gas away from the disc under the con- 
ditions valid in many galactic haloes. In the absence of a 
virial shock, the gas is not heated to the virial temperature 
until it falls all the way to the disc, where the collapse stops 
and the gas is heated in a thin layer. This may alter some 
of the assumed processes of disc formation and in particular 
the star formation rate in it. It may work against blowout 
by supernova-driven winds in dwarf galaxies. The result of 



heating near the disc instead of at the virial radius may re- 
sult in weakening the soft x-ray emission from such haloes 
and producing a high flux of La instead. In this paper we 
evaluate the conditions for the existence of a virial shock in 
galactic haloes. 

Initial density perturbations are assumed to grow by 
gravitational instability, reach maximum expansion, and col- 
lapse into virial equilibrium at roughly half the maximum- 
expansion radius. During the initial phase, and roughly until 
shells start crossing each other near the virial radius, the gas 
pressure is negligible compared to the gravitational force, so 
the shells of gas and dark matter move in a similar manner. 
Once interior to the virial radius, where shells tend to cross 
and the gas density becomes high enough, the gas pressure 
becomes an important player in the dynamics. Its hydrody- 
namic properties allow transfer of bulk kinetic energy into 
internal energy and the pressure prevents gas element from 
passing through other gas elements and from being com- 
pressed without limit. This makes the infall velocity vanish 
at the centre. Since in the cold infalling gas the typical veloc- 
ity is higher than the speed of sound, the information about 
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this inner boundary condition cannot propagate outwards in 
time, and these supersonic conditions create a shock. After 
the gas crosses the shock, it is heated up, the speed of sound 
increases, and the flow becomes subsonic. 

The shock transfers the kinetic energy that has been 
built during the collapse into internal gas energy just behind 
the shock. A stable spherical shock would slowly propagate 
outwards through the infalling gas, leaving behind it hot, 
high-entropy gas that is almost at rest. The temperature of 
the post-shock gas roughly equals the virial temperature. 
The persistence of the shock depends on sufficient pressure 
by the post-shock gas, which supports it against being swept 
inwards due to the gravitational pull together with the in- 
falling matter. Radiative gas cooling makes the gas lose en- 
tropy and pressure, which weakens the pressure support be- 
hind the shock front. Our approach here is to evaluate the 
existence of a virial shock by analyzing the gravitational sta- 
bility of the supporting gas behind the shock in the presence 
of significant cooling. 

In §2 we first summarize the standard analysis of an 
adiabatic shock and then generalize the gravitational stabil- 
ity criterion to the case where cooling is important. In §3 
we describe our spherical hydrodynamic Lagrangian code, 
which includes gravitating dark-matter and gas shells, arti- 
ficial viscosity, radiative cooling and centrifugal forces. We 
test the code in this section and in Appendix A. In §4 we 
apply the numerical code to simulations which demonstrate 
the shock formation and test the validity of the analytical 
model. In §5 we apply the shock stability criterion to re- 
alistic haloes forming in cosmological conditions. In §6 we 
summarize our results and discuss potential astrophysical 
implications. 



2 SHOCK STABILITY ANALYSIS 

Our goal here is to derive a criterion for the existence of a 
virial shock in terms of the properties of the infalling gas 
just in front of the shock front. It is based on a gravitational 
stability analysis of the post-shock gas. We first remind our- 
selves of the standard stability analysis in the simple adia- 
batic case, and then derive a more general criterion for sta- 
bility in the radiative case, under certain assumptions and 
using a perturbation analysis. 



2.1 The standard adiabatic case 

Throughout this paper, we treat the baryons as an ideal 
monoatomic gas. Their equation of state could therefore be 

written as 



(7 - l)ep. 



(1) 



where P is the pressure, e is the specific internal energy, p 
is the density of the gas and 7 is the adiabatic index. Along 

an isentrope (an adiabatic process of constant entropy) the 
pressure and density are related via P cx , so the adiabatic 
index is defined by 



d\nP 
d\np 



For a monoatomic gas 7 = 5/3.^ 

The virial shock is assumed to be a spherical accretion 
shock which propagates outwards slowly while infalling gas 
crosses it inwards. 

The kinetic energy of the infalling gas is transformed at 
the shock front into thermal energy — the post-shock gas is 
thus heated to a temperature close to the virial temperature 



of the system of dark-matter halo and gas, V-^ 



Because the original temperature of the infalling gas is neg- 
ligible compared to the virial temperature, the system obeys 
the strong-shock limit. When we denote the pre-shock and 
post-shock quantities by subscripts and 1 respectively, the 
jump conditions across the shock are in this case (Zel'dovich 
& Raiser 1966): 
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where u stands for radial velocity, Ms is the shock velocity, 
Na is Avogadro's number, Na/ p is the average number of 
molecules per unit mass, and fc_g is Boltzman's constant. 

According to standard shock theory, the post-shock gas 
is always sub-sonic (in the frame of reference of the moving 
shock) because of the increase of the sound velocity behind 
the shock. This gas is thus capable of providing the necessary 
pressure to support the shock against the gravitational pull 
inwards applied by the self-gravity of the gas and the dark- 
matter halo as well as the pressure applied by the infalling 
matter at the shock front. 

The criterion for gravitational stability of this post- 
shock gas in the adiabatic case is the standard Jeans star 
bility criterion: 7 > 4/3 (e.g., Cox 1980, Chapter 8) 

If the post-shock gas is gravitationally unstable, it falls 
into the galaxy centre on a dynamical tiincscalc and can no 
longer support the shock. As a result, the shock weakens 
and it is swept inwards. 

The Jeans criterion can be qualitatively understood 
in terms of the following heuristic derivation. For a shell 
of radius r, we compare the gravitational pull inwards, 
ag = GM/r^ (where M is the mass interior to r), to the 
pressure pushing outwards, Op = p~^VP. We assume an 
isentrope, P cc p'^ . We also assume homology, such that the 
local density scales like the mean density in the sphere in- 
terior to r, p oc M/r^ . Then VP can be replaced by ~ P/r 
and we obtain 



oc p 



7-4/3 



(7) 



(2) 



^ As the temperature exceeds the binding energy of the hydrogen 
and helium atoms, electrons become detached from tiic nuclei 
and 7 becomes smaller. Once the gas becomes fully ionized, the 
original value of 5/3 is restored, but with a different effective 
density. This should have only a marginal effect on our results, 
and is ignored in this paper. 
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If 7 < 4/3, we have an unstable configuration. Starting in 
hydrostatic equilibrium, Op/ag = 1, a perturbation involv- 
ing contraction is associated with a larger p, and therefore 
flp/flg < I by oq. (7), implying that the pressure cannot pre- 
vent collapse. If 7 > 4/3, the pressure force increases until it 
balances the increased gravitational pull. We note that even 
this simple derivation of the Jeans criterion had to assume 
homology — an assumption that we will have to adopt also 
in our analysis of the radiative case below. 



2.2 Shock stability under radiative cooling 

We wish to replace the adiabatic Jeans criterion by a more 
general stability condition that will be valid also in the rar 
diative case. This criterion must depend on the cooling rate 
and should therefore be naturally oxprossod in terms of time 
derivatives. We generalize the adiabatic 7 of eq. (2) by an 
effective 7 following a comoving volume element along its 
Lagrangian path: 



7eff = 



dlnP/dt 
d\n.p/dt 



(8) 



We expect that the system would be stable when -yes is 
larger than a certain critical value, the analog of the re- 
quirement 7 > 4/3 in the adiabatic case. 

In our Lagrangian analysis all the quantities (r, u, p, P, 
etc.) refer to comoving shells; they are all functions of the 
gas mass m interior to radius r and time t. Derivatives with 
respect to time following a comoving volume clement will be 
denoted by an upper dot, and derivatives with respect to m 
will be denoted by a prime. 

The effective gamma can be related to its adiabatic ana- 
log given the cooling rate and other post-shock gas quanti- 
ties. The time derivative of eq. (1) yields: 



P={f-l)iep + ep). 



(9) 



Energy conservation in the presence of radiative losses can 
be expressed by 
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where q is the radiative cooling rate [to be discussed below, 
e.g., eq. (20)] and V = p~^ is the specific volume. Substi- 
tuting e from eq. (10) into eq. (9), and using it in eq. (8), 
we obtain 
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Note that in the limit q/e <C p/p we reproduce the adiar 
batic case; the process is nearly adiabatic when the cooling 
timescalc is long compared to the contraction timoscalc. 

We assume that in the region close behind the shock the 
pattern of the velocity field is homologous. By this we mean 
that at any given time the (radial) velocity is proportional 
to the radius (as in a Hubble flow) , namely 



u/r = ui/ts , 



(12) 



thus providing a boundary condition for the post-shock gas. 
The homology is shown to be a valid approximation in the 
simulations discussed below, where the post-shock shell tra- 
jectories axe roughly parallel to each other in the log r — t 
plane, at any given time close enough to shock crossing. 



The time evolution of the density can then be evaluated via 
the continuity equation in Lagrangian form for the spheri- 
symmetric case. 
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where the last equality results from the assumed homology, 
eq. (12). The homology thus implies that p/p at a given t is 
a constant in m throughout the post-shock region. Eq. (11) 
can then be simplified: 
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We start with a hypothetical unperturbed state for the 
post-shock gas, where we assume that the net force vanishes, 
r = 0. The system adjusts itself to this state on a timescale 
associated with the speed of sound Cs, provided that it is 
much higher than the infall velocity u. This is expected to 
be the case in the sub-sonic post-shock medium, where Cs 
becomes high and u becomes low. The unperturbed equation 
of motion in Lagrangian form is then 

r = -47rrV- ^ =0, (15) 

where M is the total mass interior to radius r. 

We then introduce a perturbation due to a homologous 
infall velocity u. Over a short time interval 5t, it introduces 
a small displacement inwards, 5r = u 5t. In order to distin- 
guish between stability and instability wc wish to determine 
whether the induced acceleration, 5r, is positive or nega- 
tive, tending to decrease or increase the velocity respectively. 
Note that under homology, eq. (12), the relative displace- 
ment is 

6r/r = ui St/r-s . (16) 

Writing the equation of motion, eq. (15), but for the 
perturbed quantities P + 5P and r + 5r, and subtracting the 
unperturbed eq. (15), we obtain to first order 

4GM 5r 



Sr = -4nr'^{6Py + 



(17) 



Wo next manipulate the right-hand side of eq. (17) to obtain 
a simple expression involving jcS- 

In the second term we use the homology, eq. (16), and 
then the unperturbed equation of motion, eq. (15), to obtain 
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The manipulation of the first term is somewhat more 
elaborate. Wo use the definition of 7eff, eq. (8), to write 



5P = (p/p)P'y,s5t = -{3ui/r,)P-y,sSt, 



(19) 



where the second equality is due to eq. (13). Note that the 
m dependence in this term is only in the product P'yeS- 
We now express 7off in terms of the cooling rate q as in 
eq. (14), and need to take the derivative (Pq/e)' . We make 
here the standard assumption that the radiative cooling rate 
is proportional to density, 



g = pA{T) , 



(20) 



with A(T) the macroscopic cooling function and T the post- 
shock temperature. The immediate post-shock medium is 
assumed to be isothermal, reflecting via the jump conditions 
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an assumed approximate uniformity of the pre-shock gas 
over a short time interval. Using eq. (1) we have P/e = 
(7 — l)p, and together with eq. (20) it becomes 



Pq/e = (7 - l)Ap2 



(21) 



In the computation of (Pq/e)', we first replace p' = 
{dp/dP)P' , then use eq. (8) to write dp/dP = p/P, use 
eq. (21) backwards to replace (7 — l)Ap^/P by q/e, and fi- 
nally use eq. (14) to obtain {Pq/e)' = 3Mi/rs[— 27^"'(7 — 
7eff)P']. We thus have in the first term of the rhs of eq. (17) 



SuiSt 



P'b - 27cfr (7 - 7off)] ■ 



(22) 



With the right-hand side of eq. (17) given by eq. (22) 
and eq. (18), the first-order equation finally becomes, 
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Since ui and P' are both always negative, the desired sign 
of Sr is determined by the sign of the expression inside the 

square brackets. Note that in the adiabatic case, g = 0, wc 
have 7eff ~ 7, so wc recover the standard stability criterion, 
7 > 4/3. In the radiative case, 7eff ^ 7, we finally obtain 
the generalized stability criterion: 

27 



7eff > 



7 + 2/3 



— O'crit 



(24) 



For a monoatomic gas, where the adiabatic value is 7 = 5/3, 
the threshold for stability is 7crit = 10/7 = 1.43, which is 
close but not identical to the adiabatic threshold 4/3. 



2.3 Stability in terms of pre-shock quantities 

Next, we wish to express 7eff and the stability criterion in 
terms of the properties of the pre-shock gas; the infall ve- 
locity Wo and the gas density po at Vs. We use the jump 
conditions, eq. (3) through eq. (6), in eq. (14). In eq. (4) 
we assume Ug = 0, namely that the shock is temporarily at 
rest, which should be valid when the shock is marginally sta- 
ble (or unstable). This is because a stable shock is pushed 
outwards by the post-shock gas, while cooling reduces the 
pressure, slows the outward motion, and eventually causes 
it to halt and then be swept inwards by the infalling mat- 
ter and gravitational pull. The transition from stability to 
instability can thus be associated with a transition from ex- 
pansion to contraction of the shocked volume. 
According to eq. (1) and eq. (5) we have 
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According to eq. (13) and eq. (4) with Us = we have 
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With these and eq. (20) we obtain the desired expression for 
the effective 7 of the post-shock gas in terms of the pre-shock 
conditions: 



7eff = 7 - 



(7-^1)" por.A(ri) 



6(7-1)2 |moP ■ 
For a monoatomic gas, 7 = 5/3, we obtain 



(27) 



7eff = 5/3 - 18.96 
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(28) 



Based on eq. (24), the criterion for stability of a 7 = 5/3 gas 
finally becomes 



< 0.0126 . 



(29) 



The post-shock temperature is related to the pre-shock 
infall velocity using the jump condition, eq. (6), which for 

7 = 5/3 gives 



Ti 
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2 
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(30) 



For a given cooling function A(r), eq. (29) is a simple cri- 
terion for determining whether a stable shock can form at 
some radius rs of the halo. It is in a form that can be directly 
tested against hydrodynamic simulations (§3), and can serve 
for evaluating shock stability under realistic conditions in 
cosmological haloes (§5). 

Under the simplifying assumption that the gas is un- 
clumped, the cooling rate is given by eq. (20). The macro- 
scopic cooling function A(T) is related to the microscopic 
Amic(2^), the energy- loss rate of a particle, via A(T) = 
{NaX^ / n'^)^mic{T), where x is the number of electrons per 
particle. Wc assume a Helium atomic fraction of 0.1 for //, 
and X- The microscopic cooling function is shown in Fig. 1 
for three different values of mean metallicity Z. The cool- 
ing at temperatures below lO^'K is very slow because the 
main available cooling agent is molecular hydrogen, which 
is very incfliciont. At temperatures slightly above lO^K the 
cooling function peaks duo to La emission from atomic hy- 
drogen. At very low mctallicities, a second peak arises near 
lO^K due to recombination of atomic Helium. Metals give 
rise to a higher peak at ~ lO^K and slightly above, due to 
lino omission from the heavier atoms. At ~ lO^K and above, 
the cooling is dominated by brclimstralung, and the cool- 
ing function increases slowly. We use the cooling function 
as derived by Sutherland & Dopita (1993), and presented in 
their table, in the manner described in Somerville & Primack 
(1999). 



3 THE SPHERICAL HYDRO CODE 

We test the validity of the shock stability criterion using 
numerical simulations based on a spherical hydrodynamics 
code which follows the evolution of shells of dark matter and 
gas. Since the problem we intend to examine is of global 
spherical symmetry, and since we need to follow the cooling 
and the shock with high precision, we use a one-dimensional 
code. Most of the simulations presented here were run using 
2000 gas shells and 10, 000 dark-matter shells. A compara- 
ble resolution in a three-dimensional code would require on 
the order of 10^^° and lO'^^ particles respectively, which is 
impractical. We use no smoothing in the dark-matter shell- 
crossing scheme, wo introduce small-scale smoothing at the 
halo centre to avoid an artificial singularity there, and we 
include small artificial viscosity in the hydrodynamics. Tests 
of the code performance are described in Appendix A. 
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Figure 1. The microscopic cooling function of low density gas 
for three different mean metallicities as indicated in solar units. 
The peaks are dominated by atomic H and He (for Z = 0) and 
by line emission from heavier atoms (for the higher Z values) 



3.1 Dark matter 

The dark-matter particles are represented by infinitely thin 
spherical shells of constant mass and of radii r that vary 
in time. The shell of current radius r obeys the equation of 
motion 



a r 



(r- 



+ 



(31) 



where M and m refer respectively to the mass of dark mat- 
ter and gas within the sphere of radius r. The last term is a 
centrifugal acceleration, determined by the the specific an- 
gular momentum j of the particle represented by the shell. 
This j is assigned to each shell at the initial conditions and is 
assumed to be preserved during the simulation. The param- 
eter a is the smoothing length that becomes effective only 
near the centre; it has been set to be 50pc throughout this 
work. The dark-matter shells are allowed to cross each other 
(and the gas shells). The dark-matter mass is evaluated by 



M(r) = ^ AM, + i ^ S(r = r,)AM, , 



(32) 



where the shell radii and masses are denoted by rt and AA/i 
respectively, i = 1, rid- The second term adds half the 
mass of a shell when r coincides with one of the shells. Gen- 
erally, this summation requires calculations for the dark 
matter alone. The particles are kept sorted by radius. When 
two shells cross each other, we re-sort the array by exchang- 
ing pairs which violate the order. This kind of sorting algo- 
rithm, termed 'Shell's Method' in Numerical Recipes (Press 
1997), is natural in cases where only a few shells cross each 
other in each timestep. When two shells cross, they exchange 
an energy of G AAIi AMj jr. In order to conserve energy, the 
radius at which the shells cross must be known with great 
precision. We therefore reduce the timestep to a small value, 
fsc, when two shells are about to cross each other (see be- 
low). 



3.2 Gas 

The hydrodynamic part of the code is based on Lagrangian 
finite elements in the form of spherical shells. The basic 
equations governing the dynamics of each shell are 



G{M + m] 
{r + af 



dm 



P= (7-l)ep, 



(33) 
(34) 

(35) 
(36) 



An artificial viscosity term, a, is added to the pressure 
for numerical purposes, as explained below. The smoothing 
length effective at the center, a, is the same as for the dark 
matter, eq. (31). As in the model described before, the loss 
of internal energy due to radiative cooling is represented by 
the cooling rate q. 

The gas is divided into discrete shells. The mass en- 
closed within a shell. Am, is assumed constant in time, while 
the inner and outer shell boundaries move independently in 
time. Each boundary is characterized by a temporal position 
r, velocity v and specific angular momentum j. The acceler- 
ation [eq. (33)] is evaluated at the boundary position. The 
variables p, P, e, g, and T for each shell are evaluated within 
the shell between the boundaries. 

In particular, the pressure term in eq. (33) is evaluated 
at the outer boundary of shell i using eq. (35): 



--V{P + a) 
P 



Am 



[{P + a),+i ~ {P + a),] 



(37) 



where Am = (Am^ -I- Ami+i)/2. 

The boundary conditions for the outer boundary of the 
system are P — a = 0, and zero mass beyond the outer 
boundary. 

Since gas shells cannot cross each other, the gas mass in 
the sphere interior to each gas shell is constant throughout 
the simulation: 



i(rO = ^Ar 



(38) 



For the evolution of the dark-matter shell at r, we evaluate 
the gas mass that appears in eq. (31) using 

i-l 



i{r) = J2 + 



„3 



-Ami 



(39) 



i=i 



where i refers to the gas shell for which r^-i ^ r < ri 



3.3 Integration and Timestep 

The discrete integration of r and v is performed by a Runge 
Kutta fourth-order scheme (Press 1997). The state of the 
system at the beginning of each timestep is kept in memory 
until the timestep is completed, such that it is possible to 
return to the beginning of the timestep and retry with a 
smaller timestep if the convergence criteria are not met. The 
timesteps are set such that the position r and velocity v 
do not change by too much during a single timestep. For a 
given accuracy parameter trk , we demand that the difference 
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between the forth-order displacement Ar4 and the analogous 
first-order displacement An obeys |Ar4 — Ari|/r < erk, both 
for the dark-matter and the gas. The similar requirement 
is applied to the change in velocity over a timestep. If this 
condition is not fulfilled, we reduce the timestep by a certain 
factor and repeat the calculation over this timestep. We use 
here as our default £rk = 0.1. 

In addition, we make sure the timestep for each shell 
does not violate the Courant condition, for an accuracy 
parameter tc- This implies CsAf/Ar < Cc, where = 
{dP/dp)s = ^P/p is the speed of sound. We use here as 
our default ec =0.3. 

A third limitation on the timestep comes from the desire 
to conserve energy when shells cross. When two shells are 
about to cross each other within the current timestep dt, we 
set the timestep to min(dt,isc), and keep it small until they 
actually cross. We use here as our default tsc = 10-*Gyr. 

The values for ec,erk and tsc were chosen empirically 
such that energy is conserved and the dynamics converges 
to our satisfaction, in the sense that it does not change by 
much when smaller parameters are used. We demonstrate in 
Appendix A how well these requirements are met. 

Once we have computed the new radii and velocities of 
the shells at the end of the timestep, we correct the energy 
of the gas for the —PdV work term using the states of the 
system at the beginning and at the end of the timestep. The 
cooling is explicitly subtracted from the internal energy af- 
ter the hydrodynamic timestep is completed. Once the final 
state of the system is ready, it is copied onto the memory 
array of the initial state, and the simulation is ready to ex- 
ecute a new timestep. 

3.4 Initial conditions 

The simulation starts at high redshift, z = 100, with a small 
spherical density perturbation. The initial density fluctusr 
tion profile is set to be proportional to the linear correlation 
function of the assumed cosmological model, representing 
the typical perturbation under the assumption that the ran- 
dom fluctuation field is Gaussian (see Dekel 1981, and Ap- 
pendix C). The amplitude of the density fluctuation at the 
initial time, averaged over a given mass, determines the time 
of collapse, as desired. The initial velocity field is assumed to 
follow a quiet Hubble flow and the radial peculiar velocities 
build up in time. We assume the standard ACDM cosmology 
with Qm = 0.3, = 0.7, h = 0.7 and 0-8 = 1. 

3.5 Angular momentum 

We assume that in a real system the orbits of dark-matter 
particles, and the initial orbits of the gas particles, are quite 
elongated. Cosmological N-body simulations show that the 
velocity distribution tends to be more radial than tangen- 
tial (Ghigna et al. 1998; Safran & Dekel 2003) and already 
for an isotropic distribution the eccentricities are about 1:6. 
The processes we study in this paper occur away from the 
galactic disc at a radius on the order of the virial radius, 
namely in a regime where the centrifugal force can be ex- 
pected to be negligible compared to the gravitational force 
and the gas pressure force. The prescribed angular momen- 
tum for the shells is thus mainly for numerical purposes, to 



avoid divergent densities of gas or dark matter shells when 
they pass through the halo centre. Our results concerning 
the virial shock are insensitive to the actual way by which 
we assign angular momentum to each shells. 

In the current study we practically assume that the 
dark-matter particles are almost on radial orbits. The angu- 
lar momentum of the gas is prescribed such that the shells, 
once they lost their energy by radiation, would settle into an 
exponential disk with pure circular motions and a character- 
istic radius of a few kpc, smaller than the inner characteristic 
radius of the halo. Our spherical 'disc' thus contains gas that 
is cold and dense compared to the shocked gas. 

3.6 Artificial viscosity 

It is impossible to follow the discontinuous behavior across 
the shock using the conventional continuity equation for 
the density and standard conservation of energy and mo- 
mentum. The jump conditions can be calculated explicitly, 
as in eq. (3) to (6) (termed 'the characteristic method' or 
'Godonov's method'). Alternatively, as proposed by Von- 
Newman, one can slightly smear the discontinuities and then 
solve them within the framework of the standard hydrody- 
namic equations. By adding an artificial pressure term in a 
few shells around the shock, the differential equations be- 
come solvable and one can continue the calculation without 
affecting the energy and the dynamics of the shock (while 
its internal structure naturally changes). 

Artificial viscosity is applied when the inner and outer 
shell boundaries at ri and r2 approach each other, Av = 
V2 — vi < 0, and when the volume of the shell decreases, 
dV /dt = 47r(r2i'2 — rivi) < 0. The artificial viscosity then 
takes the form 

a = a2p{/^vf -\- aipCs\Av\ . (40) 

The quadratic, common form of artificial viscosity smears 
discontinuities over about 3 shells. The Linear discontinuity 
affects a slightly larger range, and is usually added with a 
smaller coefficient ai. The coefficients ai and a2 are varied 
for different shells in the course of the simulation in order 
to overcome a specific rmmerical problem in the cold 'disc', 
where the gravitational and centrifugal forces balance each 
other and the pressure force is negligible. In this case the 
gas is not a standard hydrodynamic gas because the pres- 
sure does not regulate large discontinuities, and information 
is not transported because of the low sound speed. When 
a 'disc' shell vibrates, it is artificially heated by the artifi- 
cial viscosity in every contraction until its pressure grows 
and stops the process. If wo arc not careful to properly tunc 
the artificial viscosity we may end up with one 'disc' shell 
that has been heated to lO^K while the rest of the 'disc' is 
at lO^'K. This imposes an undesired drastic decrease in the 
corresponding timestep. In order to overcome this numeri- 
cal problem, we gradually turn off the quadratic term (02) 
of the artificial viscosity inside the 'disc'. We define a 'disc' 
radius -Rdisc to be the largest radius for which the difference 
between the gravitational and centrifugal forces is less than 
1/4 of the gravitational force. Once at r < O.QRdmc, we con- 
tirmously decrease the parameter 02 in cq. (40) according to 
02 = (r — 0.3i?disc)/(0.3-Rdisc) and make it completely vanish 
at r < 0.3i?disc- This prescription was found by trial-and- 
error to properly solve the numerical problem in most cases. 
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Figure 2. Simulation of the adiabatic case. The sequence of (solid 
red) curves describe the log radii of Lagrangian gas shells as a 
function of time. The simulation was of 2000 gaseous shells (shown 
here) and 10, 000 dark matter shells. The radius of every 20th 
gaseous shell is plotted. Shown on top are the virial radius and 
the shock. The shock exists at all times. It gradually propagates 
outwards, and it practically coincides with the virial radius. 



The linear term of the viscosity, being proportional to the 
speed of sound, is anyway very small in the cold 'disc', so 
effectively no artificial viscosity is applied in the inner 'disc'. 

Appendix A provides tests and examples of the hydro- 
dynamic simulations in some detail. 



4 VIRIAL SHOCK IN THE SIMULATIONS 
4.1 Existence of a virial shock 

We now investigate the formation of a virial shock using the 
spherical hydrodynamical simulations described above. We 
wish to test in particular the validity of the analytic stability 
criterion developed in §2. 

In order to mimic a typical perturbation in a ran- 
dom Gaussian field (Dekel 1981, Appendix B), the initial 
density-fluctuation profile was set to be proportional to the 
correlation function, normalized such that the mean den- 
sity fluctuation in a sphere enclosing Mi — 10^^ Mq was 
Si = 0.09 at z — 100. For example, the shell encompassing 
M ~ 3 X 10^" Mq is expected to collapse at z = 3, and 
M ~ 10^^ Mq is expected to collapse at z — 0. 

Fig. 2 shows the time evolution of the radii of La- 
grangian shells in a simulation of the adiabatic case, with 
the cooling turned off. We find that a shock exists at all 
times. It appears as a sharp break in the fiow lines, associ- 
ated with a discontinuous decrease in infall velocity [eq. (4)]. 
Shown in the figure is the shock radius, defined by the outer- 
most shell for which the inner and outer shell boundaries ap- 
proach each other and the volume of the shell decreases [the 
same conditions that have been used for turning on the arti- 
ficial viscosity in eq. (40)]. The shock gradually propagates 
outwards, encompassing more gas mass and dark matter in 
time. The gas below the shock is pressure supported and 
at quasi-static equilibrium. Not shown here are the dark- 
matter shells, which collapse, oscillate and tend to increase 
the gravitational attraction exerted on the gas shells. 



Shown in comparison is the evolution of the virial ra- 
dius, computed from the simulation density as the radius 
within which the mean overdensity is Av times the mean 
cosmological background density. The virial overdensity Av 
is provided by the dissipationless spherical top-hat collapse 
model; it is a function of the cosmological model, and it 
may vary with time. For the Einstein-deSitter cosmology, 
the familiar value is Av — 176 at all times. ^ For the family 
of flat cosmologies {flm + ~ 1), the value of Av can be 
approximated by Bryan & Norman (1998) 



Av ~ (IStt^ + 82a; - 39a;^)/(l 



(41) 



where x = Qmiz) — 1, and ^Imiz) is the ratio of mean matter 
density to critical density at redshift z. For example, in the 
ACDM cosmological model that serves as the basis for our 
analysis in this paper (ilm = 0.3, Qa = 0.7), the value at 
2: = is Av — 340. We see in Fig. 2 that the shock radius 
almost coincides with the virial radius at all times. This is 
hardly surprising, as the shock is likely to appear at the 
outermost radius at which shell crossing first occurs, which 
is near the virial radius (to be demonstrated in Fig. 8 below) . 

Fig. 3 is the result of a similar simulation, but now with 
realistic radiative cooling for Z = 0. We see that a stable 
shock does not exist in this case before t = 3.9Gyr. During 
this period, the cooling makes the gas lose its pressure sup- 
port and lets it collapse freely under gravity into the halo 
centre. The collapse stops by the assumed angular momen- 
tum, in a 'disc' whose marked radius can be identified at 
the bottom of the plot by the abrupt change of the infalling 
flow lines into horizontal lines. The matter in the 'disc' is 
angular-momentum supported. As is visible in the figure, a 
shock, in the sense of a discontinuity in velocity and density, 
is present at the edge of the disc. Once the stability criterion 
is met, a shock forms and propagates outward abruptly. The 
propagation of the shock causes it to re-enter a regime for 
which 7cfi is below the critical value. Consequently, the post 
shock gas becomes non-supportive again and falls. The os- 
cillatory behavior of the shock continues with an increasing 
period until it stabilizes at the largest radius for which the 
stability criterion is met. The shock never expands beyond 
the virial radius because shells do not tend to cross there 
(Fig. 8 below). 

Fig. 4 shows the evolution of the total mass interior to 
the characteristic radii in the adiabatic and radiative simu- 
lations of Fig. 2 and Fig. 3. In the adiabatic case, the shock 
mass practically coincides with the virial mass, and the gas 
never forms a disc. In the radiative case, the shock is ini- 
tially at the disc radius, and only after 3.9Gyr does it start 
propagating outwards. 



4.2 Testing the stability criterion 

We first test the assumption made in §2.2 that the post- 
shock gas has a homologous velocity profile in the vicinity 
of the shock. In Figures 2 and 3 we see that the flow lines 
beneath the shock tend to be nearly parallel in the log r — t 



^ This can be derived from the top-hat formalism of Appendix B, 
once the final radius is assumed to be fixed at half the maximum- 
expansion radius but the overdensity is evaluated at the time 
when the top-hat sphere would have collapsed to a singularity. 
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Figure 3. Simulation of the radiative cooling case, with Z = 0. 
The curves are as in Fig. 2, with the 'disc' radius added. There 
is no shock outside the 'disc' at early times, when the virial mass 
is small, because the cooling is too efficient. A shock develops at 
later times, when the mass is larger, and it quickly propagates out- 
wards. After a couple of oscillations the shock radius approaches 
the virial radius. 
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Figure 4. Evolution of the total mass interior to the virial radius, 
shock radius and 'disc' radius, in the adiabatic simulation (top) 
and the radiative simulation (bottom). 



Figure 5. Model predictions for shock stability. Top: The flow 
lines describe a simulation similar to that shown in Fig. 3, ex- 
cept that the cooling rate has been set to be unrealistically 
high such that the shock cannot develop. The 'disc' radius is 
marked by the break in the flow lines at the bottom. Shown 
on top are 4 contours of equal 7(,ff values (dashed lines, for 
7;,ff = —1.0, 1.0, 1.428, 1.55) They are computed by the model, 
[eq. (27)], from the local ( "pre-shock" ) gas quantities correspond- 
ing to the flow lines at the background, but assuming realis- 
tic cooling in the computation of 7cff . The long-dashed (green) 
curve corresponds to 7eff = 7crit — 1.428, above which a shock 
would have formed under realistic cooling. The negative value, 
7off = —1, is not physical. Middle: Time evolution of 7eff at infall 
just above the 'disc' radius (solid) in comparison with the model 
critical value for stability 7crit = 1.428. Since 7(,fl is monotoni- 
cally increasing with decreasing radius, the shock is expected to 
form first at the 'disc' radius. Shock formation is predicted by 
the model at t ~ 3.9Gyr. Bottom: Evolution of the characteris- 
tic masses in the simulation with the unrealistically high cooling 
rate. The virial mass at the time of shock formation is predicted 
to be about lO^^M©. 



plane. This means that dlogr/dt = u/r ~const., namely 
homology, eq. (12). A similar behavior has been found in all 
the simulations that we have performed. 

Next, we use the simulations to test the validity of the 
stability criterion derived in §2, eq. (29), where jcti is given 
by eq. (27). Recall that the 7eff is expressed in terms of the 
pre-shock quantities. In order to map the value of 7off in 
the different regions of the free falling gas in the forming 
halo, we ran the same simulation as in the previous section 
except that the coohng rate was set to be very high, such 
that the virial shock never develops. Fig. 5, top panel shows 
the flow lines in this case, on which overlaid are 4 contours 
of equal 7off values, evaluated via eq. (27) with eq. (30) and 
the cooling rate from Sutherland & Dopita (1993). As shells 



are falling into the halo, their 7off is gradually increasing. 
Also, as time progresses, the value of 7eff at the same radius 
is increasing. By following the value of 7eff just above the 
"disc" radius (shown as the break at the bottom of the plot), 
in comparison with the critical value for stability 7crit, we 
can therefore use our model to predict when we expect the 
virial shock to form. This is shown in the middle panel of 
Fig. 5. We see that at early times (and smaller masses) we 
have 7cff < 7crit, predicting no stable shock. The system 
is predicted to enter the stable-shock regime at about t — 
3.9Gyr, where 7cff becomes larger than 7crit. A comparison 
with the realistic radiative simulation described in Fig. 3 
yields that this model prediction is very accurate: the shock 
indeed starts forming at i = 3.9Gyr, and is globally stable 
thereafter. 
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Figure 6. Time evolution of ■y^g at the 'disc' radius or the shock 
radius, whichever is larger (top), and the associated characteristic 
masses (bottom), in the simulation with realistic cooling of Fig. 3. 
The shock forms for the first time when 7^3 becomes larger than 
7crit = 10/7. The shock then oscillates while it approaches a 
steady state near the virial radius. 



Fig. 6 shows the actual evolution of 7off at either the 
'disc' radius or the shock radius, whichever is larger, as com- 
puted directly from the pre-shock quantities in the simula- 
tion with realistic cooling. Shown at the same times are the 
characteristic masses, already shown in Fig. 4. The virial 
shock is first generated at time t = 3.9Gyr, when the virial 
mass is 10^^ Mq. Starting at this time, the shock is propagat- 
ing outwards very rapidly. As a result of this fast expansion, 
the 7cff of the pre-shock infalling matter at the shock, which 
is decreasing with r, drops below the threshold. This makes 
the shock lose its pressure support, it becomes temporarily 
unstable and its expansion slows down until it is eventually 
swept back on a dynamic time scale. The associated drop in 
total mass behind the shock, seen around t = 5.8Gyr, is due 
to the fact that the dark matter is not swept back with the 
gas. Once the shock is shrunk to a low enough radius, ^^s 
rises again to above 7crit; the shock becomes stable again 
and it resumes its associated expansion towards the virial 
radius. After the conditions for the shock stability are first 
met, the shock is visible most of the time. In the rest of 
this paper, we treat haloes at this state as ones containing 
a virial shock. 

Fig. 7 presents results similar to Fig. 6 for two other 
simulations with different initial overdensities, and therefore 
different masses collapsing at different times. The small dif- 
ference seen in one case between the 7eff at which the shock 
actually forms and the predicted 7crit may be totally due 
to numerical inaccuracies in the simulation. Such inaccura- 
cies may occur when a dark-matter shell crosses a gaseous 
shell, which, near the threshold, may lead to a slightly pre- 
mature shock formation. This is seen in the convergence test 
of our code described in Appendix A, when we compare the 
shock formation times in table A2. Thus, the model stability 
criterion is found to be valid within the accuracy of the sim- 
ulations in all the cases studied, indicating that the model is 
not limited to a special range of masses and collapse times. 
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Figure 7. Same as Fig. 6 for two other initial conditions. In 
all 3 cases the cosmology is ACDM and the initial fluctuation 
density profile is proportional to the correlation function, but 
with different amplitudes, as indicated a.t z = 100 within the 
sphere enclosing 10^^ Mq. Note the different time scale in the 
bottom panel. The qualitative behavior of the shock is similar in 
the three cases, and so is the success of the model in predicting 
shock stability. 



5 SHOCK STABILITY IN COSMOLOGY 

The analysis of §2 thus provides a successful criterion for 
shock stability, eq. (29), as a function of the pre-shock prop- 
erties of the infalling gas at radius r: the density, velocity 
and metallicity. In order to apply this criterion to a given 
protogalaxy in a cosmological background, we wish to eval- 
uate the gas density and velocity just before it hits the disc, 
for a gas shell initially encompassing a total mass M that 
virializes at redshift Zv In this calculation we assume an 
Einstein-deSitter cosmology, as a sensible approximation at 
z > 2 (where Q.^ > 0.9). 

We assume a given universal baryonic fraction /b , and a 
global spin parameter A which determines the ratio of disc to 
virial radius. The initial mean density perturbation profile, 
Si{M), is given at some fiducial time in the linear regime; 
it is the average profile derived from the power spectrum of 
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initial density fluctuations, as described in Appendix C. In 
the cosmological toy model used here we approximate the 
power spectrum as a power law, Pk oc k" , where n ~ —2.4 
to mimic the ACDM power spectrum on galactic scales. 

Wc follow gas shells from the initial perturbation till 
they approach the disc using a two-stage model. During the 
expansion, turn-around, and until an assumed virialization 
at half the maximum-expansion radius, we assume no shell 
crossing, the total mass interior to the shell remains con- 
stant in time, and wc follow the radii, density and velocity 
of the shell via the spherical top-hat model (see Appendix 
B). Prom the virial radius inwards we assume that the gas 
shells, which do not cross each other, contract inside the 
fixed potential well of an isothermal dark-matter halo. This 
idealized model involves several crude approximations, such 
as the instantaneous transition at the virial radius, and ne- 
glecting the effect of the angular momenta of the individual 
gas particles at small radii, but we show using spherical sim- 
ulations that this model predicts the minimum halo mass 
for which a stable shock first appears to an accuracy bet- 
ter than 25%. This allows us to use the model for exploring 
the critical mass as a function of cosmological parameters 
such as galaxy formation time, metallicity, spin parameter, 
fluctuation power spectrum, and baryonic fraction. 



5.1 Toy model until virialization 

For a given shell M and initial mean perturbation profile 
5{M) [standing for the <5i(M) of Appendix B], the top-hat 
model [eq. (B8) and eq. (B2)] yields the implicit solution 



r{M,r)) = rv(M) (1 - cos?;) , 
t(M,77) =U{M) (77 -sin 77), 



(42) 
(43) 



where the mass dependence enters via the virial quantities 
= CiM'/^<5(M)-\ (44) 



U = rv/vv = Cf/^G"'/^ 5{M) 



-3/2 



(45) 



with = GM/r^. The coefficient d = (6/7r)^''^(0.15/pui) 
is determined by pui, the cosmological density at the initial 
time when 5{M) is given, independent of M. 
The velocity of the shell M is 



dr/dri 
dt/dr) 



sm?7 



(1 — cos rf) 



(46) 



At virialization, r] = 3tv/2, it is simply u = —Vv 

In order to evaluate the local density, we follow the radii 
of two adjacent shells, encompassing masses M and M + 
dM respectively, at a given time t, e.g., the time when shell 
M virializes (at half its maximum expansion radius). Let 
?7 correspond to shell M at that time, and ri + drj io shell 
M -h dM. In order to express dr) in terms of dM we use the 
fact that the time t is the same for the two shells: Q = dt = 
{dt/dM)dM + [dt/drfjdr]. Using eq. (43) this gives 



3 (n — sin 77) 5' , 
2(1 — cos jy) 



(47) 



where we denote 5' = dd/dM. Expressing dr in terms of dM 
and dr) based on eq. (42), we obtain using eq. (47) and after 
some algebra 



dr _ 1 dM / _ 3Md' 
r ~ 3 M \ S 



1 - 



3 sinr;(77 — sin 77) 

2 (l-COS7?)2 



(48) 



At virialization of shell M, rj = 3it/2, the quantity in square 
brackets equals (10 + 9n)/4. Not surprisingly, if the initial 
perturbation is of uniform density, 5' = 0, we are left with 
dr/r = {l/3){dM/M), the straightforward result of M cx r^. 
Recall that the virial radius rv of shell M can be obtained 
either from the universal density at the time of virializa- 
tion using eq. (B4), or from the initial perturbation using 
eq. (B8). The desired local density p can be obtained from 
eq. (48) via dM = Aur'^pdr. 

If the initial perturbation profile is a power law, 5{r) cx 
using M DC we have 3M5' /5 = -{n + 3). So 

finally 



(49) 



eq. (48) [or cq. (49) in the power-law case] allows us to 
compute the desired radii of the two adjacent shells at the 
time of virialization of shell M. 



5.2 Toy model after virialization 

Given the radii of the two adjacent shells at rv, we enter 
the shell-crossing regime and continue to follow the shells 
down to the disc radius by numerical integration. The shell 

radius r and velocity u are related via energy conservation. 
Assuming that the gas shells contract without crossing each 
other inside a dark-matter halo that is a fixed isothermal 
sphere, the total mass interior to the shell that originally 
encompassed a total mass M is 



M(r) = fi,M + {l- ,M—r, 

and the gravitational potential at r is 



0(r) = -^-(l-/.)^ln(-) 
r Tv \ r / 



GM . 



(50) 



(51) 



The integration is performed by advancing r according 
to the velocity u and then recalculating u according to en- 
ergy conservation: 

(l/2)w2 + tp(r) = {l/2)v^ + <P{r^) = const. (52) 

We follow shell M for the time it falls from r = t\- to the 
disc radius r = Afv, and shell M + dM for the same time 
interval. Denoting the separation between the shells at the 
end of this time interval by dr, we compute the desired gas 
density by 

/b dM 

47rr2 dr 

The resultant values of r, u and p are inserted into eq. (28) 
in order to obtain an approximation for 7eff and then to 
evaluate stability by eq. (24) . This allows us to check stabil- 
ity for the case where mass M virializes at redshift Zy, with 
metallicity Z, spin parameter A, baryonic fraction /b, and a 
given power spectrum. 

5.3 Model versus simulations 

In Fig. 8 and Fig. 9 we compare the evolution of the quanti- 
ties of a given gas shell according to the toy model described 




Figure 8. Model versus simulation. Evolution of the pre-shock 
gas properties along a trajectory of one Lagrangian shell as a 
function of the radius. The solid (red) curve corresponds to the 
smallest shell from the simulation of Fig. 3 for which the shock 
formed. The dashed (green) curve is the prediction of the toy 
model, calibrated here in an ideal way to match the simulation 
at maximum expansion. The bottom panel shows the -y^ff derived 
from the above quantities for the model and the simulation, in 
comparison with the critical value of ■y^ff marked by the horizontal 
line. We notice that indeed -y^fi = ■jait as r ^ 0. The fit between 
model and simulation is remarkable. 



in the previous subsections and according to the spherical 
hydro simulation described in the earlier sections. We fol- 
low a specific shell that hits the disc at about t ~ 3.8Gyr, 
just before the shock starts propagating into the halo [see 
Fig. 3]. The quantities shown as a function of radius r are 
total mass M interior to r, radial velocity u, gas density p, 
and the corresponding value of 7cff . For M and u the evo- 
lution starts at the top-left corner and ends at the bottom 
left, while for p the upper part of the curve corresponds to 
the expansion phase and the lower part to the contraction 
phase. The evolution of 7 is followed only during part of the 
contraction phase. 

In Fig. 8 we calibrate the toy model to match the simu- 
lation at the maximum expansion radius. We see that while 
the mass interior to the shell is reproduced by the model 
only to a limited accuracy in the last stages of the collapse, 
the velocity, density and the resulting value of 7eff are re- 
covered very well by the model. This allows us to predict 
quite accurately the point where 7eff exceeds 7crit. 

Since we wish to use the toy model without an exact 
knowledge of the conditions at maximum expansion, we nor- 
malize the model evolution in Fig. 9 based on M and Zv 

The slight deviations in u and in p now translate into 



Figure 9. Model versus simulation. Same as Fig. 8, except that 
the model was calibrated in a practical way based on M and Zv 
without reference to the simulation results. The model predic- 
tions now deviate somewhat from the simulation results, but the 
deviation amounts to less than 25% in the estimate of the critical 
mass. 

a larger error in 7cff . The error in the toy model originates 
mostly from the slight ambiguity in the definition of the 
virial radius. On one hand we assume it to equal half the 
maximum-expansion radius, and on the other hand we as- 
sume it to represent an overdensity of Av as in eq. (41). 
These two assumptions are not fully consistent with the ac- 
tual behavior of the virializing system in the simulation. 
Nevertheless, we see below that our approximate model al- 
lows us to estimate the critical halo mass below which the 
shock does not form to an accuracy of better than 25%, 
which is quite satisfactory for our purpose here. 

5.4 Critical mass for shock formation 

Fig. 10 shows for several different cases the critical halo 
mass, below which a shock does not propagate into the halo, 
versus the redshift at which this critical mass virializes. For 
each case we compare the model prediction to the shock 
formation as actually seen in the simulation. The cases dif- 
fer by the mean metallicity, Z — and Z — 0.05 for the 
lower and upper sets of points respectively, and by the am- 
plitude of the initial perturbation, corresponding to a range 
of shock-formation redshifts at every given Z. The assumed 
baryonic fraction is always /b — 0.13, but the assumed spin 
parameter may be different for the different shells in a given 
simulation because we set it for each shell such that the fi- 
nal disc has an exponential surface density profile. However, 
the A values vary in the range 0.02 to 0.05, compatible with 
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Figure 10. Model (green circles) versus simulations (red 
squares). Displayed is the critical halo mass for shock formation 
versus the virialization redshift. In all cases = 0.13, n = —2.4, 
and 0.02 < A < 0.05. The high set of points is for Z = 0.05 and 
the low set is for Z = 0. 



the distribution of spin parameter in cosmological simula- 
tions (Bullock et al. 2001). We see that the model predicts 
the critical mass with an accuracy better than 25%, such 
that we can use it for mapping the parameter space in more 
detail. 

Fig. 11 shows for several different choices of parame- 
ters the critical halo mass for shock formation versus the 
halo virialization redshift as predicted by the model. A virial 
shock does not form in haloes of masses below the line. The 
lines are not always monotonic due to the non-monotonic 
features in the cooling curves (Fig. 1). Shown in compari- 
son is M*, the characteristic mass for haloes forming at z 
according to the Press-Schechter approximation (Lacey & 
Cole 1993). The default values of the parameters, used un- 
less specified otherwise, are Z = 0, \ = 0.05, /b = 0.13 and 
n = -2.4. 

The upper panel has the metallicity varying from Z — 
to Z — 0.3. The critical mass tends to be higher at higher 
redshifts (especially for Z = 0) because the higher density 
implies more efficient cooling. It is striking that even for the 
case of zero metallicity, for which the cooling is not at its 
maximum efficiency, an halo cannot produce a shock 
until a relative late redshift, z ~ 2.1. The addition of a 
small amount of metals, Z = 0.05, increases the cooling rate 
significantly (see Fig. 1) such that haloes start producing 
virial shocks only after z ~ 1.6. 

The second panel has A varying as marked. The shock 
forms slightly earlier if the disc is smaller (lower A), because 
the conditions become more favorable for shock formation 
closer to the centre. At high redshifts the increase in infall 
velocity happens to balance out the increase in density, tem- 
perature and cooling rate. The post shock temperature there 
is a few lO'^K. The bottom line is that the critical mass is 
not too sensitive to A. 

The third panel has /b varying as marked. The critical 
mass is monotonic with the baryonic fraction because the 
cooling rate is monotonic with gas density. The parameter 
/b can be interpreted as the fraction of the baryons that ac- 
tually take part in the shock formation. This can be smaller 
than the universal baryonic fraction if some of the gas falls 
into the halo in the form of dense clumps. Even with /b as 




Figure 11. Critical mass for shock formation as a function of 
virialization redshift, according to the model tested against sim- 
ulations in Fig. 10 (dashed curves). Shown for comparison is the 
characteristic press-Schechter M* as a function of z for ACDM 
(solid curve). The default values of the parameters are: /b = 0.13, 
A = 0.05, Z = 0, and n = —2.4. One parameter is varied in each 
panel as indicated. 



low as 0.05, meaning that most of the gas is not participat- 
ing in the cooling, an halo would not produce a shock 
until z ~ 2.4. The conclusion is that the critical mass is not 
too sensitive to /b either. 

The bottom panel explores three values for the initial 
power index n approximating the power spectrum of ACDM 
on galactic scales. The dependence of the critical mass on n 
in this regime is weak. 



6 DISCUSSION 

The heating of the gas behind a virial shock in haloes has 
been a basic component in galaxy formation theory (Rees & 
Ostriker 1977). We studied the conditions for the existence 
of such a virial shock in spherical haloes. We first pursued 
an analytic stability analysis in the presence of cooling, and 
then demonstrated its validity using high-resolution spheri- 
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cal hydrodynamical simulations. The obtained criterion for 
shock stabiUty in terms of the post-shock quantities is 



7eft- = /, ,^ > 10/7 



(54) 



din p/dt 

In terms of the pre-shock gas properties, this condition reads 
PorsA(Ti) 



wo ■ 



< 0.0126 , 



(55) 



where po and mq are the gas density and infall velocity in 
front of the shock, rs is the shock radius, A(T') is the cooling 
function which depends on the metallicity Z, and Ti oc Mq 
is the post-shock temperature as a function of the pre-shock 
infall velocity. 

Based on this criterion, we find that a virial shock forms 
only in big haloes forming at late redshifts. A virial shock 
does not form in smaller haloes forming early where the cool- 
ing behind the shock efficiently removes its pressure support. 
For example, we find that most galactic haloes that have 
collapsed and virialized by a: ~ 2 did not produce a virial 
shock. Haloes loss massive than ~ 10" M© never produce a 
shock even if the gas is of zero metallicity. If the metallicity 
is non-negligible (e.g. Z ~ 0.05), this lower bound to shock 
formation rises to ~ 7 • 10^^ Mq. When a shock does not 
exist, the gas is not heated to the halo virial temperature 
until it falls all the way to the disc at the inner halo. 

Forcada-Miro & White (1997), in an unpublished work, 
have pursued independently a numerical analysis along sim- 
ilar lines, involving a more detailed treatment of the cooling 
processes involved. They also find that the virial shock ra- 
dius is significantly reduced due to the cooling in haloes of 
small masses, M < 10" Mq. In their case the shock never 
completely disappears because of a different feature in their 
numerical scheme; they put all the cooled post-shock gas in 
one central "shell" to avoid numerical difficulties at the cen- 
tre. This makes the inner boundary of the system follow the 
shock quite closely in cases where there is efficient cooling 
behind the shock, and allows the presence of a small-radius 
shock even in such cases. Overall, our numerical results are 
in encouraging agreement, and our analytic model provides 
a natural explanation for their numerical results as well. 

The most severe uncertainty when attempting to apply 
our results to real galaxies arises from the assumed spherical 
symmetry in both the model and the numerical simulations. 
The validity of this approximation for the asymmetric halo 
configurations in the hierarchical clustering scenario is an 
open question to be addressed in future work. Nevertheless, 
we notice that Katz et al. (2001) and Fardal et al. (2001) 
observe in their cosmological simulations that a large frac- 
tion of the mass accreted onto haloes indeed remains cold 
and is never heated to the virial temperature. Toft et al. 
(2002) find in their simulations, using a similar treeSPH code 
as the one used by Katz et al. (2001), that the soft X-ray 
radiation is mainly emitted from within the innermost 20kpc 
of their haloes, well inside the virial radius, in encouraging 
agreement with our results. On the other hand, it is not ob- 
vious that the resolution in these simulations is adequate for 
studying the shock physics involved; our estimates indicate 
that three-dimensional simulations with proper resolution 
are not practical at present (Appendix A). 

Another complication may axise from radiative effects. 
Even when there is no virial shock, the kinetic energy of 



the gas eventually turns into radiation when the gas in- 
fall motion is brought to a halt at the disc. At such den- 
sities, the width of the shock front is much smaller than the 
width of the cooling front behind the shock, ~ 10^'^pc ver- 
sus ^ lO^pc. Thus, the gas in a thin shell behind the shock 
is heated to a temperature corresponding to its kinetic en- 
ergy, and it cools by radiating soft X-rays. The X-ray rar 
diation is expected to generate an ionized Hn bubble, in 
which the ionization rate balances the recombination rate. 
The Stromgren radius of this bubble is relatively small, on 
the order of a few kiloparsecs, because the high gas density 
implies a high recombination rate. The recombination pro- 
cess then generates a flux of La radiation, emitted at the 
inner few kpc of the halo. 

A naive inspection of cross sections might indicate that 
the La radiation would be trapped inside the halo. This 
could in principle affect the shock stability in three difi'erent 
ways: by increasing the radiation pressure, by heating up the 
infalling matter, and by slowing down the radiative cooling 
responsible for the shock instability. It has been argued by 
Rees & Ostrikor (1977) that the radiation pressure at these 
low temperatures must be insignificant compared to the gas 
pressure even if all the internal energy was drained from the 
baryons into the radiation field. One might add that since 
the radiation pressure behaves like ay — 4/3 gas, it could at 
most make the system marginally stable. When work is done 
on the radiation field, any leakage of radiation out would 
turn the energy into cooling rather than PdV, and will thus 
reduce the effective gamma, making the system unstable. 

Partial heating of the infalling gas should not affect our 
analysis as long as the temperature of the infalling matter is 
significantly below the virial (post-shock) temperature such 
that the strong shock approximation remains valid. The ef- 
fect of the reduced cooling rate is yet to be investigated. In 
practice, we do not expect the radiation trapping to be very 
efficient, because the effective opacity is reduced by ther- 
mal broadening and by the systematic blue shift due to the 
gas infall motion. When the opacity is high, the radiation 
heats up the gas, which enhances the thermal broadening 
and the coUisional ionization rate. This reduces the opac- 
ity and allows for radiation escape. The system is likely to 
reach a steady state in which it gradually cools. This process 
is under current investigation. 

Feedback effects may further complicate the picture and 
affect shock stability. The energy fed back to the gas from 
stars, supernovae and AGNs may heat the halo gas and ex- 
pel part of it. Merging substructures may have additional 
complicated effects. These effects cannot be captured by our 
idealized spherical analysis, and a proper study would prob- 
ably require high-resolution three-dimensional hydrodynam- 
ical simulations. While observations and certain theoretical 
considerations indicate that feedback effects are likely to be 
important in galaxies as large as ~ 10^^ M0 and may thus 
affect the shock stability (e.g. Dekel & Silk 1986; Dekel & 
Woo 2003, and references therein), it has proven difficult 
for numerical simulations to reproduce such effects in any 
but small dwarf galaxies, indicating that feedback effects 
may not be so crucial for the understanding of shock sta- 
bility. Until the dust settles on the role of feedback effects, 
our preliminary conclusions based on the spherical analysis 
should be taken with a grain of salt. 
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The general absence of a virial shock might have three 
direct implications, which we study in associated papers. 

First, as explained above, when the g£is is heated at the 
disc rather than near the halo virial radius, the generated 
X-ray radiation serves to ionize the gas and is not emitted 
outwards. The result would be a suppression of the X-ray 
emission in the range 5 x 10^ to 2 x lO^K. This may help 
explaining the missing X-ray problem pointed out by Pen 
(1999) and Benson et al. (2000). Pen (1999) argue that there 
is an ordor-of-magnitude discrepancy between the soft X-ray 
flux as observed by Cui et al. (1996), after subtracting the 
contribution of quasars, and the predicted flux from haloes 
constructed by a Press-Schechter hierarchical model under 
the assumption of shock heating to the virial temperature. 

Second, the infall energy, via the ionizing X-ray, is effi- 
ciently transformed into La radiation at the inner few kpc of 
the halo. A related increase in the La flux has indeed been 
seen in the cosmological simulations of Fardal et al. (2001). 
This may explain the observed high flux of La emitters at 
high redshift (e.g. Pentericci et al. 2000, 2001; Breuck et al. 
2000) . Based on the high observed flux and the assumption 
that the La is emitted from stars, Pentericci et al. (2000) 
estimate large masses for the La emitters, but the much 
higher flux per unit mass predicted by our model may lead 
to significantly lower mass estimates. Based on our analysis, 
most of the La flux is expected to be emitted from the inner 
few kpc of the halo, where the gas is at ~ 10'*K. Neglect- 
ing line shifts and broadening, the halo might be opaque to 
La , thus eventually emitting its energy from an outer pho- 
tosphere where the halo becomes transparent. However, a 
careful study of the thermal broadening and the systematic 
redshifts within the halo is required in order to determine 
whether the system is opaque or transparent to the La pho- 
tons. This is a subject of an ongoing investigation. 

Finally, the direct collapse of cold gas into the disc may 
have interesting theoretical consequences to be worked out. 
It may induce an efficient star burst in analogy to the burst 
originating in the shock between two colliding gas clouds. In 
turn, the strong inwards flow of gas may prevent an efficient 
gas removal by supernova-driven winds. In particular, cur- 
rent cosmological semi-analytic models (SAM) of galaxy for- 
mation (Kauffmann et al. 1993; Kaffmann et al. 1999; Cole 
et al. 1994; SomervUle & Primack 1999; Mailer et al. 2001, 
and related works) use the standard picture of heating be- 
hind a virial shock in their modeling. This has strong effects 
on the disk formation rate, star formation rate, feedback 
etc. Other semi-analytic models (Efstathiou 2000; White & 
Frenk 1991) also appeal to the slow gas infall rate as a mech- 
anism that regulates the gas input into the disc. Since the 
cooling time for a 10^^ Mq halo is relatively short, the SAM 
predictions for such haloes may be only slightly affected by 
the inhibition of heating. However, given some metal en- 
richment, no heating is expected for haloes as massive as 
~ 7 X 10^^ Mq, for which the cooling time is longer, and the 
effect on the SAM predictions may be more severe. Shocks, 
when present, are also expected to alter the properties of 
the gas, for example - extinct dust particles. These effects 
can change SAMs that incorporate dust extinction. 
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APPENDIX A: TESTING THE HYDRO CODE 

The numerical code, Hydra, has been developed specifi- 
cally for simulating the evolution of a single spherical halo 
through collapse and feedback processes. A proper compu- 
tation of the cooling and shock formation requires high pre- 
cision. In this appendix we describe a few of the tests per- 
formed in order to veriiy that the code works properly. In 
the following three subsections we test for energy conserva- 
tion, spatial convergence, and the performance of the code 
in a self-similar case. 



Al Energy Conservation 

Our numerical scheme does not use the total energy equa- 
tion in the integration of the partial differential equations. 
Furthermore, the total energy of the system is not a straight- 
forward sum of other variables that are involved in the cal- 
culation. The requirement of energy conservation is there- 
fore an independent test for the accuracy of the numerical 
scheme. Energy conservation is harder to achieve than spa- 
tial convergence for several reasons. First, the error in total 
energy is systematic, in the sense that when dark-matter 
shells cross each other the energy tends to increase. Second, 
since our system is only marginally bound, the total energy 
is a small difference between two large quantities. We notice 
that energy conservation is simpler to achieve when there is 
no cooling, or when dark matter is absent (and thus there 
is no shell crossing). 

The total energy of the system at time t is the sum of 
terms: 



Table Al. Energy conservation 



E = Ka+Ta + Kg + Tg + U + Q, 



(Al) 



where subscripts d and g refer to dark matter and gas re- 
spectively, K stands for kinetic energy, T stands for poten- 
tial energy, U is the gas internal energy, and Q is the thermal 
energy lost to radiation by time t. 

For the dark matter, these are straightforward sums 
over the discrete dark-matter shells: 



Ki = AMi 



i=l ^ 

^ ^GAMM^ 



2 _|_ Ji 



(A2) 



(A3) 



where A4i is the total mass interior to dark-matter shell i, 
as defined in §3. 

For the gas shells, recall that the quantities r, v and j 
are given at the inner and outer shell boundaries, i — 1 and i 



description rig 




tsc/Gyr Bfin/Binit 


nominal 


2k 


10k 0.3 0.1 


10-* 0.992 


small e's 


2k 


10k 0.1 0.01 


10-5 0.994 


more shells 3k 


15k 0.3 0.1 


lO"'' 0.997 


Table A2. Shock formation times and energy conservation 


ngas 




^ final / ^initial 


formation of shock(Gyr) 


3,000 


15, 000 


0.997 


3.9 


2,000 


10, 000 


0.994 


3.9 


1,000 


5,000 


0.976 


4.01 


500 


2,500 


0.923 


3.76 


250 


1,250 


0.711 


3.27 


125 


625 


0.252 


2.82 



respectively, so we compute the shell energies by averaging 
over the two boundary values: 



ViT-f + Vi-iri_i 



Kg= - ^ A nn 

1 "g / • 2 , . 2 \ 2 

2 ^ V + '"i-i 
^ [(r? + r3_,)/2]i/3+a 



(A4) 



(A5) 



The internal energy is a straightforward sum 
U = ^ a Ami . (A6) 
The energy radiated away, Q = J dt J qdm, is computed by 



(A7) 



j=i i=i 



where At-' is the length of timestep j, and ql is the cooling 
rate in shell i at timestep j (in units of erg g^^s^^). 

In a run with 10,000 dark-matter shells and 2,000 gas 
shells, we require and obtain energy conservation at the 
level of 1% in a Hubble-time, using a typical Runge Kutta 
timestep of about 5 x IQ-^Gyr. (Such a run takes about 10 
hours on an Alpha-6 DEC processor). 

We check the conservation first by varying the accu- 
racy parameters presented in §3.3, and then by varying the 
number of shells. The three cases presented in Table Al 
demonstrate that the results converge when the accuracy is 
increased. The simulations shown in this table are of the 
standard case with realistic cooling shown in Fig. 3. When 
cooling is shut off, energy conservation is much better. With 
the nominal choice of accuracy parameters the final energy 
is 0.9999 of the initial energy. 

A2 Spatial Convergence: 3D versus ID 

A proper treatment of the competition between the pres- 
sure increase due to contraction and the pressure decrease 
due to cooling requires high temporal and spatial resolution. 
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Figure Al. Self-similar profiles by our hydro code with gas only. 
Top: Density profiles at times 4 to llGyr. Bottom: The same 
profiles normalized according to the self-similarity law r ac t^^^ 
to t = 5Gyr. 



In particular, when the spatial resolution is increased, the 
shock appears earlier. Table A2 shows results from simula- 
tions of the case with realistic cooling (Fig. 3), all with the 
same accuracy criteria (ec, trk, and tsc), but with different 
spatial resolutions. The average distance between gas shells 
near the center ranges from about 80pc to 2kpc. With the 
poorest resolution of 125 gas shells the virial shock appears 
almost immediately after the virialization of the first shells 
of the simulation. The energy changed by about 75% during 
this simulation. Even if we assume that the precision of a 3D 
calculation is as good as that of an analogous ID calcula- 
tion (actually SPH codes converge slower than finite element 
schemes for problems involving shocks), we still need to cube 
the number of particles or grid points in order to achieve the 
same resolution. A three-dimensional simulation with 2 x 10® 
gas particles and 2.5 x 10* dark-matter particles, which is 
close to the limit of what is computationally feasible today, 
would correspond to the unsatisfactory case with the lowest 
spatial resolution in table A2. 



A3 A self-similar case 

When the initial conditions are scale free (unlike the initial 
conditions assumed in the body of this paper, motivated 
by ACDM) , and when the cooling function is also scale free 
(unlike the realistic cooling function used above), the results 
should be self similar. This can provide a test for the accu- 
racy of our numerical code. We follow Bertschinger (1985a,b) 
in using an initial perturbation consisting of a point-mass 
embedded in a uniform-density background. Far from the 
point mass, the system should be self similar. We ran a sim- 
ulation of such a case using our code with gas only (/b = 1) 
and no cooling, starting at 2 = 200 with an overdensity of 
10 inside the innermost 2kpc. 

The upper panel of Fig. Al shows the density profile at 
different times. As expected, a shock appears at every time 
as a density jump by a factor of 4 [= (7 -I- l)/(7 — 1) for 7 = 
5/3], and the post-shock gas settles to a complete rest after 
it is shocked. (The slope of the post-shock density profile 



-26 

-28 

-30 
-26 

-28 

-30 



0.&Gyt^:r 


Original profiles 


0.9Gyr ""^"-..^^^^ 




IGyr 




1 1 n,y,r 

1 . 1 Oyi 




2Gyr 






Normalized to 1 .3Gyr 







10 



100 



radius [kpc] 



Figure A2. Similar to Fig. Al but with dark matter included, 
/b = 0.13. The profiles are scales to t = l.SGyr. 



is somewhat different from Bertschinger (1985b), because 
our calculation assumes a ACDM cosmology rather than the 
Einstein-deSitter assumed by Bertschinger.) The lower panel 
shows the same profiles after they were scaled to the same 
time (5Gyr) according to the scaling relation of Bertschinger 
(1985b): r <x t^^^. We see that our simulations recover the 
expected scaling relation almost perfectly. 

Fig. A2 shows an analogous test for the case where both 
gas and dark matter are present, with /b = 0.13. The results 
are similar except for the somewhat higher noise level caused 
by the dark-matter component. 



APPENDIX B: TOP-HAT MODEL 

Consider a bound spherical perturbation encompassing mass 
M, whose mean density fluctuation profile at some fiducial 
initial time in the linear regime is Si{A4), embedded in an 
Einstein-deSitter (EdS) cosmological background when the 
universal expansion factor is ai . We wish to express the shell 
radius r as a function of time in terms of M. 

The implicit solution for a closed "mini-universe", via 
a conformal time parameter rj specific to this perturbation, 
is 

r = rv(l — COS77) , (Bl) 
t^—i-n-sm-n). (B2) 

Vv 

Maximum expansion occurs when r; = tt, and then a collapse 
to half this radius, which we identify with virialization, is ob- 
tained at 77 = 3n/2, with virial radius and corresponding 
virial velocity v'^ — GM /r^. We normalize the universal ex- 
pansion factor a by identifying it at the initial time with the 
shell radius r. Assuming rn <^ 1, we have n ~ (l/2)rv?7? and 
ti ~ (l/6)(rv/«v)77f , so a; = n yields a; = {Qr^vl /2Y^^t1'^ . 
The EdS expansion factor, a oc t^^^ , can now be related 
using eq. (B2) to the perturbation's r) at any time: 



= {9r,vl/2f'h^''-' = (9/2)i''Vv(r? - sinr^f^^ . 



(B3) 



The mean density within the perturbation relative to the 
universal density at the same time becomes a straightfor- 
ward function of rj: 
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_ 9 (r; - sin 77)^ 



Pu 2(1 — cos 77)3 



(B4) 



This is a standard result of the top-hat model. 

In order to relate the density to the small initial per- 
turbation at rji -C 1, we obtain from eq. (B4) by a proper 
Taylor expansion to the first non-vanishing order: 



5i ~ 0.15??i , 



(B5) 



where 6 = p/pu — 1 is the mean fluctuation. Using this in 
the linear term of eq. (Bl) we obtain 

Oi = (l/2)rv??? = (l/0.3)rv^i . (B6) 

This allows us to write the mean density at any ij, using 

eq. (Bl), as 



a; 



Sf 



Pui r3 0.3(1 - cos 7?)3 ■ ^^"^^ 
Recalling that p = M/[(47r/3)r^] we finally obtain at any rj 

Ml/3 



r = Ci 



(1 - cosr?) , 



where 



Ci 



/6y/3 ai5 

y-rrJ Pui 



(B8) 



(B9) 



In particular, at 77 = 3n/2, we obtain for the virial radius 
rv = CiM^^^/5i{M). The constant d is independent of M; 
the universal density pui [= (1 + 2:i)^Puo] is determined by 
the choice of the fiducial redshift Zi at which 6i{M) is given. 



APPENDIX C: INITIAL PROFILE 

We adopt in the linear regime the typical density fluctuation 
profile for the assumed power spectrum of fluctuations. For 

a Gaussian random field, this profile is proportional to the 
two-point correlation function (Dekel 1981): 



5iW = <5oi|(oy, 



(CI) 



where Soi specifics the amplitude normalization. For a given 
power spectrum P{k), the correlation function is given by 
(Peebles 1993, eq. 21.40) 



= 47T / 

Jo 



kr 



(C2) 



and the local variance is 



^(0) = 47r / k^dkP{k) . (C3) 
Jo 

The mean density fluctuation interior to radius r, con- 
taining mass M = (47r/3)pui?'^ when the fluctuation is small, 
is 



Mr] 



(C4) 



This involves the integral (Peebles 1993, eq. 21.62) 



Mr) ^ 



f 

Jo 



3 Jo 



k'^dkP{k)Ws{kr), (C5) 



where Ws{kr) is the Fourier transform of the top-hat win- 
dow. 



W,{kr) = 3 



sin kr cos kr 



(kr)^ (kr) 



(C6) 



with ji the spherical Bessel function. 

We use in the simulations of this paper the ACDM 
power spectrum as from the fitting formula of Bardeen et 
al. (1986) 

AkT'^ik), 



P{k) 



(C7) 
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T{k) = {1 + [ak/r + {hk/Tf'^ + (cA;/r)Y} 

with a = 6.4/i-iMpc, b = S-Oh'^Mpc, c = l.T/i-^Mpc, 
u = 1.13 and F = 0.21 (the rCDM model of Efstathiou et 
al. 1992). The normalization is such that as ~ i- 

In the cosmological toy model we approximate the 
ACDM by a power- law power spectrum Pk oc fe", for 
which the two-point correlation function is also a power law, 
^(r) oc r-^"+^\ and then 
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where n provides the normalization. In terms of mass we 
obtain 



Si{M) 



M 



(47r/3)puir3 



-(n-|-3)/3 



(CIO) 



This serves as the input to eq. (B8), or eq. (48). 

We normalize the initial perturbation such that a spe- 
cific mass Ml reaches virialization at some cosmological 
epoch Ov = 1/(1 + Zv). Using eq. (B3) and eq. (B6) we 
obtain the linear analog to the nonlinear fluctuation growth 
rate: 



1/3 



At virialization, this gives 6v 

\ CLv J 



sinjj)^''^ . 
; l{r) = 37r/2) ; 



(Cll) 
1.58. Then: 
(C12) 



The normalization parameter 5oi (or n ) at Oi is obtained by 
equating this with eq. (C4) [or eq. (CIO)] at M = Mi. 



